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Abstract 

In this paper we present a high-order kernel method for numerically solving diffusion and reaction- 
diffusion partial differential equations (PDEs) on smooth, closed surfaces embedded in R d . For two- 
dimensional surfaces embedded in R 3 , these types of problems have received growing interest in biology, 
chemistry, and computer graphics to model such things as diffusion of chemicals on biological cells 
or membranes, pattern formations in biology, nonlinear chemical oscillators in excitable media, and 
texture mappings. Our kernel method is based on radial basis functions (RBFs) and uses a semi- 
discrete approach (or the method-of-lines) in which the surface derivative operators that appear in the 
PDEs are approximated using collocation. The method only requires nodes at "scattered" locations 
on the surface and the corresponding normal vectors to the surface. Additionally, it does not rely on 
any surface-based metrics and avoids any intrinsic coordinate systems, and thus does not suffer from 
any coordinate distortions or singularities. We provide error estimates for the kernel-based approximate 
surface derivative operators and numerically study the accuracy and stability of the method. Applications 
to different non-linear systems of PDEs that arise in biology and chemistry are also presented. 

Key words: radial basis functions • mesh-free • manifold ■ collocation ■ method-of-lines • pattern for- 
mation ■ Turing patterns ■ spiral waves. 
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1 Introduction 



Kernel methods such as those based on radial basis functions (RBFs) are becoming increasingly popular 
for numerically solving partial differential equations (PDEs) because they are geometrically flexible, algo- 
rithmically accessible, and can be highly accurate. Since Kansa's pioneering work [42], there have been 
many successful applications of kernel methods to various types of PDEs defined on planar regions in two 
and higher dimensions (see, for example, 20 Ch. 38-45] and the references therein). More recently, these 



methods have been developed and applied to PDEs defined on the surface of a sphere [23 25 26 30 37 
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In this paper, we present a high-order kernel method using RBFs for numerically solving certain PDEs 
defined on more general surfaces. In particular, we focus on diffusion and reaction-diffusion equations on 
smooth, closed embedded submanifolds M C M. d . In the case of two species, the prototypical form of the 
latter type of these PDEs is 

~ =S u A M u + f u (t, u, v), 

— =S v A M v + f v (t, u, v), 
at 

where u, v : M — > K, 5 U ,5 V > 0, f u , f v are (possibly non-linear) scalar functions, and Am is the Laplace- 
Bcltrami operator for the surface. For two dimensional surfaces embedded in three dimensional space, 
systems like (JlJ have received growing interest to model such things as diffusion of chemicals on biological 



cells or membranes 59 60 , pattern formations in biology 12 67 , nonlinear chemical oscillators in excitable 
media 13 ,40, 52], and texture mappings in computer graphics 65 . 

In the past decade, many methods have been developed for PDEs like on surfaces. Nearly all of these 
techniques can be classified into two types, intrinsic methods and embedded, narrow-band methods. The 
former methods use coordinates intrinsic to the surface and a surface-based mesh to discretize the differential 
operators (e.g. [Tl^[T5|fl7l[45||47||63y7l| ) , while the latter methods extend the entire PDE in R 3 in a narrow 
band around the surface and then modify the differential operators so that the solution is restricted to the 
surface (e.g. [2|[8|[9|[l4{[39 , 49 , 50 , 57 60] ) . Intrinsic methods have the benefit that the resulting discretization 
scheme is consistent with the dimension of the original problem. However, properly dealing with the inherent 
coordinate distortions or singularities that arise in the metric terms of the surface differential operators can be 
difficult, and these methods are generally limited to low orders of accuracy. Embedded, narrow-band methods 
have the benefit that the surface differential operators are posed in extrinsic coordinates so that all coordinate 
singularities can be avoided and standard methods such as finite differences on 3D Cartesian grids, or finite 
elements on 3D unstructured meshes can be used. Additionally, the surface can be naturally represented 
using standard level-set methods, and thus quite topologically complicated surfaces can be handled. However, 
these methods require consistent extensions of the initial data on the surface to the embedding space, which 
can be non-trivial. Also, some of these methods (e.g. [9j) lead to degenerate surface differential operators 
since, for example, they allow diffusion to occur only in directions tangential to the surface. This makes 
it difficult to use implicit discretization techniques for time-dependent problems. Additionally, the need to 
implement artificial boundary conditions in the embedding space can lead to accuracy degradations 39 . 
Finally, all embedded methods have an added computational expense since they solve the equations in a 
dimension at least one greater than that in which the equations are posed. This added expense can grow 
quite significantly depending on the surface, the order of the differential operators, and the order of accuracy 



of the method 49 



The present kernel method combines many of the benefits of both the intrinsic and narrow band methods. 
It uses a semi-discrete approach (or the method-of-lines) in which the surface derivative operators (e.g. 
surface Laplacian) that appear in the PDEs are approximated using collocation. These approximations 
are made using RBF interpolation on surfaces [36]. In this way, the method is similar to the ones used 
in [23,25,26,30] for the surface of a sphere. However, the new method applies to more general surfaces, 
including those that can be defined parametrically or implicitly (as the level-set of a function). Since the 
method is based on RBFs, it allows the computational nodes to be placed at "scattered" locations on the 
surface (see Figure [I] for examples). Furthermore, it does not rely on any surface- based metrics or intrinsic 
coordinate systems, thus avoiding any coordinate distortions or singularities. In this regard, our kernel 
method is similar to the embedded, narrow-band methods discussed above. The difference is that our 
method approximates the surface differential operators directly on the surface instead of having to extend 
quantities off the surface and do the approximations in K 3 . Thus, like the intrinsic methods, our method is 
posed in the same dimension as the original problem. This ability to directly approximate the differential 
operator also means that no degeneracies arise from having to modify the equations so that the solutions 
remain restricted to the surface. Furthermore, computational efficiencies are gained by not having to extend 
into the full three dimensional space to solve the problems. 

Recently, a different kernel method based on RBFs was introduced by Piret for solving PDEs on surfaces, 
which is called the Orthogonal Gradients (OGr) method [57]. This method shares many similarities to 
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the embedded, narrow band method known as the Closest Point Method 49 50 58 and is different from 
our method in a number of ways. For example, the OGr method requires expanding into the embedding 
space (]R 3 ) to obtain discrete surface differential operators, which increases the computational complexity. 
Further, a differentiation matrix is constructed by enforcing that the first and second normal derivatives 
to the surface of a certain RBF interpolant vanish. This means in particular that one must analytically 
calculate the derivatives of the normal vectors to the surface, which requires that a smooth representation of 
the surfaces is known. In contrast, our method only requires nodes at scattered locations on the surface and 
the corresponding normal vectors to the surface, and does not require expanding into R 3 . Furthermore, the 
discrete operators for the present method are constructed from the inversion of a single traditional positive- 
definite kernel interpolation matrix, which is guaranteed to be invertible. Additionally, we provide a theory 
for the convergence of the discrete surface differential operators. 

The application of kernel methods (including RBFs) to problems on general surfaces is still in its infancy. 
The goal of this paper is thus to investigate some of the theoretical and applied aspects of our kernel method 
and to establish its applicability. As a result, we (a) present error estimates for approximating various 
surface differential operators; (b) numerically investigate the stability and accuracy of the method, showing 
how solutions can be computed to very high orders of accuracy; and (c) apply the method to two relevant 
problems in biology and chemistry. 

The remainder of the paper is organized as follows. In next section, we briefly review kernel interpolation 
with RBFs. We follow this with a short presentation in Section [3] on how to formulate surface differential 
operators in Cartesian coordinates. Section[4]describes how these surface differential operators are discretized 
in the form of differentiation matrices and presents the kernel method in mcthod-of-lincs form. Making heavy 
use of the results from [36] , we provide error estimates for the discrete surface derivatives in Section [5] (and 
proofs in Appendix [A| . Section [6] numerically investigates the eigenvalue stability of the kernel method. We 
numerically demonstrate the convergence the method for the forced scalar diffusion equation to two different 
surfaces in Section [7j Applications of the method to simulations of pattern formations in a Turing system 
and to spiral waves in excitable media are presented in Section [8j We conclude with some comments on 
potential future enhancements of the method in Section [9] 

Before continuing to the main body of the paper, we pause to make some relevant remarks on what 
follows. While there are many types of classifications of surfaces or manifolds, in our work, we use the term 
surface or manifold to always refer to a smooth embedded submanifold of M. d with no boundary. We use 
the notation M to denote the manifold. For simplicity, we present our method for manifolds of dimension 
2 embedded in R 3 . However, it can be naturally extended to manifolds of codimension 1 embedded in R d . 
Although not pursued here, more technical, but straightforward extensions are possible for manifolds of 
higher codimension. 



2 Brief overview of kernel interpolation 

Fundamental to the construction of our approximate surface differential operators are kernel interpolants. 
Let f2 C M. d and 0:S!x57— s-Mbea continuous function, which we will refer to as a kernel. Given a set of 
"scattered" nodes X = {xjjj^ C fl and a continuous target function / : M. d — > M sampled at X, a kernel 
interpolant takes the general form 

N 

J /(x) = c i<H x > x i)' x e n > ( 2 ) 

where Cj are determined by requiring I$f\ x = f\ x - ^ is well-known that there are many kernels (f> for which 
a unique solution to this problem exists. In particular, if for any non-zero vector b € R N , cj> satisfies 

N N 

5^5^h i 0(x i ,x i )6 i >O 1 (3) 

then there exists a unique solution to (J2|. We call 4> that satisfy ([3| for all finite node sets X C f2 positive 
definite kernels on f2. 
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(d) Torus (e) Dupin's cyclide (f) Bretzel2 

Figure 1: Example surfaces and "scattered" nodes (represented as small solid spheres) used in the numerical 
experiments. A description of the surfaces and the node sets is given in Appendix [B| 



In this study we focus on the popular subclass of positive definite kernels called radial basis functions 
(RBFs), which have the property that 0(x,y) = 0(||x — y||), where || • || is the standard Euclidean norm. 
(Strictly speaking, there are RBFs that are only conditionally positive definite 20 Ch. 7], but our focus will 
be on the positive definite ones). For these radial kernels ^ becomes 



N 



7 /(x) = ^c^(||x-x J ||) 

3=1 

and the interpolation constraints can be expressed as the following linear system: 



(4) 







Cl 
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A x 














C N _ 







(5) 



x 



where (Ax)ij — 0(||x.j— Xj||). For a positive definite cf>, this system is positive definite and hence non-singular. 
We will use the following two popular RBFs in this study: 



Matern : <j> v {r) = C v (er) v - d / 2 K u _ d/2 (er), v > d/2, e > 0, 

IMQ: 0(r) = 1 s > 0, 
y/i + (sr) 2 



(6) 
(7) 



which are positive definite on R d . In C u = 2 1 -^- d ^/T(iy-al/2) and K p is the modified Bessel function 
of the second kind of order j3. The Matern family of kernels are piecewise smooth, with the smoothness 
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controlled by the parameter v. For example, with v = ^g^, 4> v G C 2 (E rf ) and with v = ^„ e C 4 (M d ). 
In contrast, the inverse multiquadric (IMQ) kernel is C°°(M. d ). Infinitely smooth kernels like the IMQ are 
used often in solving PDEs, and we have included results for the Matern family to illustrate the different 
convergence rates that are possible with finitely smooth kernels. 

Typically, kernel, and in particular RBF, interpolation is applied to problems where the data are scattered 
in M. d or on the surface of the sphere. The present authors have recently proved that favorable error estimates 
can be achieved when RBF interpolation is applied to reconstruction problems on more general manifolds 



M of K 36 . It is worth noting, that the interpolants for these reconstruction problems are not modified 
from their original form Q, i.e. distances are still measured as straight line distances, and not as distances 
intrinsic to the manifold. Thus, the interpolation method needs only limited knowledge of the underlying 
manifold; it only requires computing the solution to the linear system |5). In what follows, we will use Q to 
approximate the surface gradient and surface divergence over a given set of nodes X C M on the manifold, 
and combine these to obtain an approximation for the surface Laplacian. For this construction, we will only 
need a point set X C M and the normals to M at each point in X. 

We conclude with some remarks on the parameter e in (|6) and ([7]), which is called the "shape parameter" 
since it can be used for changing the kernels from peaked (large e) to flat (small e). In general, better 
accuracy is obtained for smaller values of e 20,32, Ch. 16-17]. However, the standard way of computing 
the interpolant by solving ([5) becomes increasingly ill-conditioned. While several types of stable algorithms 
have been developed for bypassing this ill-conditioning [121,27,29,31], these algorithms generally break down 
when the nodes X fall on a lower dimensional manifold of unless the surface topology is taken into 
account in the algorithm (as done on the sphere in [29]). Modifications to these algorithms are underway to 
address this issue |46 , but they are not yet generally available. As a result, in this study we will not give 
much attention to issues of the shape parameter, but will instead choose it in the applications to ensure the 
linear system ([5) is reasonably well-conditioned. 



3 Continuous surface differential operators 



Here we assume that M is a smooth embedded submanifold of K d with no boundary and dim(M) = d — 1. 
In the descriptions of the surface differential operators on M that follows, we will focus on the case of d = 3 
for notational simplicity. The extension to other d should be obvious. Additionally, all expressions are given 
in extrinsic (or Cartesian) coordinates. 

For any point x = (x,y,z) on M, we denote the normal vector to M at x as n = (n x ,n y ,n z ) and the 
vector space of tangent vectors to M at x as T X M. The surface gradient on M at x is then given as 



:= PV = (I — nn 



(8) 



Here the 3-by-3 matrix P projects vectors in M 3 to T X M. Letting e x , e y , e z , be the standard unit vectors in 
x, y, and z directions in ]R 3 , we can re- write |8) in component form as 

V] \Q X ~ 

(9) 

I 3 at a point 





~(e x 


P)V~ 




(e x — n x n) 


V" 




p x 


V" 




-gx 






P) V 






V 




p« 


V 




gy 




\e z 


P)V 




(e z — n z n) 


V 




p z 


V 




Q z 



With this notation, the surface divergence of a smooth vector field f = (f x , f v , f z ) : 
x e M can be expressed as 

v M • f := (pv) ■ f = g x f x + gyp + g z f z . 

Finally, the Laplace-Beltrami operator on M at x can be written using ^ and ([TO) as: 

A M := V M • V M = (PV) • (PV) = Q X Q X + QPQP + g z g z . 



(10) 



(11) 



In what follows, we will approximate the Laplace-Beltrami operator in the form of (11) by replacing the 
continuous operator ([9) with a discrete version based on the kernel interpolant Q. 
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4 Kernel collocation technique in method-of-lines form 



The discrete approximations to the surface differential operators that will be used in the method-of-lines 
(MOL) are based on applying the projected gradient ^ to the kernel interpolant Q. We illustrate how to 
construct these approximations by first applying the projected gradient to a radial kernel. We then show 
how to express the discrete operators as differentiation matrices. This is followed by a description of the 
MOL formulation with these discrete operators for the diffusion equation. As in the previous section, we 
formulate everything in R 3 , with an obvious extension to other dimensions. 



4.1 Projected gradient of a radial kernel 

Let x and x.,- be points on M, and let rj(x) = ||x — Xj| denote the Euclidean distance from x to Xj, i.e. 
rj(x) = \J{x — Xj) 2 + (y — Vj) 2 + (z — Zj) 2 . With this notation, a radial kernel ^ : R 3 x R 3 -> 1 centered 
at Xj is given by </>(r.,(x)). Under the assumption that (/>(r\,(x)) is differentiable at x, a simple application 
of the chain rule gives the gradient of </>(r., (x)) as 

d 



V<Xr,(x)) 



oz 



r,-(x)) 






-Xj) 




r 3 -(x)) 




(y 
{* 


-Vj) 
- z i) 


r 3 -(x)) 









(x - Xj 



r j( x ) 



(12) 



where we have used ' to denote differentiation with respect to r. Note that the apparent singularity in the 
above formula when rj(xj) = will analytically cancel since <f> is assumed to be an even function with at least 
one continuous derivative. The surface gradient of 4>(rj(x)) can then be obtained from ([9]). After applying 
([9]), the components of the surface gradient of 0(rj(x)) are 

„ ^(rj(x)) 



Q x 4 


fa(x)) = 


((1 


- n x n x )(x - 


x i) 


- n x n y (y - yj) 


- n x n z (z 


g y 4 


fa(x)) = 


((1 


— n v n v )(y — 


Vj) 


— n x n y (x — Xj) 


— n y n z (z 


g z 4 


fa(x)) = 


((1 


- n z n z )(z - 


h) 


— n x n z (x — Xj) 


— n v n z (y 



4>'{r. ] (*)) 
rj (x) 

r,-(x) 



(13) 
(14) 
(15) 



We now have all the components necessary to build the RBF approximation of the surface gradient PV of 
a scalar quantity defined on M. 



4.2 Discrete operators in terms of differentiation matrices 

Let / : M — > K be some differentiable target function known at a set of "scattered" node locations X = 
{xjl^L-L C M. The first step in constructing a discrete approximation to the surface gradient of / on X is to 
construct an interpolant to / of the form Q . In the notation of this section, the interpolant takes the form 

N 

V(x) = X;^(r J -(x)), (16) 
3=1 

where the coefficients Cj are determined by collocation. We next apply the projected gradient operator ^ 
to I<j,f and evaluate at the nodes in X. Focusing on the Q x component of this operator and using (13 1, we 
obtain 

N 

(5*ww)U = X> (^M x )))lx= v 

3=1 

- ^ c 3 ((! ~ n i n i)( x i - x j) - (»i - Vj) - n'niizi - z,)) - " < {X '" 



3=1 



r 



(Xi) 



1,3 



G 



where i — 1, . . . , N and (nf , n\, nf ) is the normal vector at Xj. We can write this approximation in terms of 
a matrix vector product using the notation from the linear system ([5]) as follows: 



(9*I+f)\ x = B x x c f = {B X A- X l )f x = G x f x . 



(17) 



The matrix G x x is an N-by-N differentiation matrix that represents the discrete RBF approximation to the 
x-component of the projected gradient operator over the set of nodes in X. We can similarly obtain the 
discrete approximations to the y- and z-components of the projected gradient operator as follows: 



{9yi^f)\ x = (B x A x ')f x = Glf x , 
(G z I <t> f)\ x = (B x A x 1 )f x = G x f x , 

where the entries in B v x and B x are given by 

( B x)i,j = (C 1 - n \ n \)^ ~ Vj) ~ n i n i(%i ~ x i) ~ n \K(. z i ~ z i)) 



(18) 
(19) 



r,-(xj) 

(Btfij = ((1 - ntnl)( Zi - Zj ) - nfnffc - xj) - nfnf ( Vi - y )) ^M*)) 



for i, j = 1, ... , N. The full discrete approximation to the surface gradient operator Vm is given in terms of 



(17H19) as the 37V-by-iV matrix 



G 



x 



G\ 
G\ 
G\ 



(20) 



A discrete approximation to the surface divergence operator (10) can also be constructed from (17)-(19|. 
Let f = (f x ,f v J z ) : M -> R 3 be a smooth vector field and f x = [/f f y x f x ] T denote the 3iV-by-l 
vector of samples of f on X. Then the discrete approximation of the surface divergence of f on X is given 

by 



[G x 



Lt x 



G' x ]t 



x 



G x x .r x 



u xJx 



Gxfx- 



(21) 



(ID 



Finally we can combine ( 20 1 and ( 21 ) to obtain a discrete approximation to the Laplace-Beltrami operator 



L x — D X G X = G X G X X + G\G\ + G Z X G Z X . 



(22) 



The N-by-N differentiation matrix L x represents our discrete RBF approximation to the surface Laplacian 
and is what we will use in the applications that follow. It is important to note that the construction of 
L x only requires a node set X C M and the normal vectors to M. Additionally, L x is constructed using 
extrinsic coordinates and is thus free from any coordinate singularities (e.g. the pole singularity for the unit 
sphere S 2 ). In the case of M = § 2 , other studies have taken a different approach at approximating the surface 



Laplacian using kernel collocation (e.g. 37 70 ). These studies apply the surface Laplacian directly instead 
of approximating it with projected gradient as we have done. For more general manifolds, it is non-trivial 
to work out the surface Laplacian in closed form, which is why we are advocating the use of the discrete 



approximation (|22 ) 
In Section 



we show that our discrete approximations (20)-(22) converge to the continuous surface 
derivatives Q-Ill), respectively, as the nodes in X increases and "fill" the manifold 



Although the computational cost of constructing the discrete surface differential operators is 0(N 3 ), these 
constructions are a preprocessing step and can be done once and stored. Applying these surface differential 
operators to a vector requires 0(N 2 ) operations. In Section [9j we make some remarks on how these cost 
may be reduced. 
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4.3 MOL formulation 



To describe the MOL formulation of our kernel method, we use the equation for diffusion of a scalar quantity 
u on a surface with a (non-linear) forcing term: 

— =SA M u + f(t,u). (23) 

Here S > is the diffusion coefficient, f(t,u) is the forcing term, and an initial value of u at time t = is 
given. Letting X = {xj}jL 1 C M and ux G denote the vector containing the samples of u at the points 



in X, our method for (23) takes the form 

d 



u x = SL x ux + f (t,u x ) , (24) 



where Lx is the differentiation matrix |22|) corresponding to our discrete approximation of Am- This is 
a system of N coupled ODEs and, provided it is stable, can be advanced in time with a suitably chosen 
time- integration method (see the Sections [7] and [8] for examples) . In Section [6j we numerically investigate 
the eigenvalue stability of the method. 



5 Convergence results for the discrete surface derivatives 

In this section we present the convergence results for our discrete differential operators. Given that these 
operators are essentially obtained by differentiating interpolants, our error estimates will rely on the ap- 
proximation power of kernel interpolants. Recently, Sobolev error estimates for a wide variety of kernel 
interpolants (constructe d fro m RBF kernels restricted to surfaces) were made available for smooth, embed- 



ded submanifolds of Mr 36 , and these results can be used to derive the bounds in this section. Precisely 
stated theorems and their proofs are given in Appendix [A] for the interested reader. A summary of the 
results will be presented after we establish the necessary notation. 

5.1 Definitions and notation 

As is typical in scattered data approximation, the error estimates are given in terms of the mesh norm h 
and separation radius q of the point set X on M, which are given by 

h := sup min d M (x, Xj) q := - min <4i(Xi,Xj), 

where cJm(x, Xj) is the geodesic distance between x and Xj on M, i.e. it is the arclength of the shortest curve 
on M connecting x and x^. The uniformity of the data is also important - this is measured by the mesh 
ratio p := h/q. 

Every positive definite kernel is the reproducing kernel for a space of continuous functions, often called 



the native space (see 20 Chapter 13] or 68 Chapter 10] for details). We will denote the native space of 
a scalar kernel <j> with domain ft C M. d by Af^Cl). When determining the approximation power of a given 
kernel, there is a rich theory to draw from when the target function lies within A/^,(f2). 

The native spaces associated with all of the kernels we consider here are subsets of, and sometimes equal 
to, Sobolev spaces. We denote by if*(R ) the space of functions whose distributional derivatives up to order 
t are all in L2(M. d ), with the usual norm. We will ultimately consider Sobolev spaces on a manifold M, where 
M is a 2-dimensional smooth embedded submanifold of M 3 . We denote by i?*(M) the space of functions on M 
whose projections onto M 2 by coordinate charts are all in the Sobolev space iJ'(K 2 ) (see (36) Section 2]). Since 
our manifold is embedded in E 3 , a natural choice for vector-valued functions on M is H'(M) := (i?*(M)) 3 . 
The norm for a function f = (f x , f v , f z ) e H*(M) is given by 

2 2 2 X 1/2 

f x \\H'(M) + WPWhuu) + II/ z IIh*(i 



H*(M) ■- \\\J \\H*(M) 
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The kernels that we use throughout the paper give rise to spaces of continuous Sobolev functions. Thus 
from this point on we assume that for some r > 3/2 the kernel <f> satisfies 



0H<C(1 + |M| 2 )- T , (25) 

where is the Fourier transform of on R 3 and C is a constant^] Radial kernels whose Fourier transforms 
decay exponentially, such as the Gaussian or IMQ, obviously satisfy this. For finitely smooth kernels, we 
will require the stronger condition that 

C (i + IMI 2 r T < £(w) < c(i + |M|T T , (26) 



where c and C are positive constants. Note that in particular the Matern kernels ^ satisfy ( 26 ) since [68| 
p. 133] 



(l + ||w 



2\ 



It can be shown that kernels satisfying (25) obey A/^(M) C H S (M) 1 and that those satisfying (26) give 



A/^(M) = iP(M) 36 Theorem 3.3], where s = r — 1/2. The parameter s governs the smoothness of the 



native spaces intrinsic to the surface M, so the error estimates that follow will often be given in terms of s. 
5.2 Summary of convergence results 

We need to define continuous analogues to our discrete differential operators. Recall that the operators Gx 
and Dx are designed by applying surface differential operators to an RBF intcrpolant. This prompts the 
following definitions: 

GmJ '■— Vm^/, Df&f ■— Vm • i$f , L M f := D m GmJ, 

where 7$f = {Isf x , I<f>f y , I<f>f z ) is the vector- valued kernel interpolant of f. For any continuous functions / 
and f , the functions above sampled on X agree with the associated discrete differential operator acting on 
the vector fx or fx, e.g. (Lm/)|x = Lxfx- To determine how well Lx approximates Am, for example, we 
will derive estimates for the error ||Lm/ — Am/|| over the entire domain. 

The error will be measured in the 2^2 (M) and L^iM) norms, and we do this for a few reasons. Even 
though error is often measured empirically with an co-type norm, the current theoretical co-norm error 
estimates for many kernel approximation problems typically give slower rates than those observed experi- 
mentally. The observed co-norm rates tend to be closer to the L2-type convergence rates, which are often 
optimal. Thus one can probably expect the faster ^(M) convergence rates in practice. With that, here is 
our first approximation result. 



Result 1: First-order discrete operators. Suppose the kernel 4> satisfies (25) with r > 3/2 + 1 and 
define s = t - 1/2. Then the following holds for all / e A^(M) and f € (N^(M))\ 



\Guf — Vm/IU 2 (M) < ), ||Gm/ - Vm/IU 00 (m) < C(h s ), 

D M f-V M -f||L 2 (M) < Oih*- 1 ), \\D M t~V M -i\\ Laam <(D(h s - 2 ). 

For a proof, see Proposition [T] Appendix |A") Note that for kernels such as the Gaussian and IMQ, which 



satisfy (25) for all r > 0, this implies convergence faster than any fixed algebraic order for target functions 
within the native space. 

A key assumption in the above result is that the target function is within A/^M). When the kernel is 
less smooth, we have the ability to provide error rates for many target functions not within A/0(M) - the 
only requirement is that the targets be in a Sobolev space smooth enough so that first-order derivatives are 
continuous (since dim(M) = 2, this means that / € H@ (M) for some /3 > 2). In this case the uniformity 
of the data sites, measured by the mesh ratio p, enters the error estimates. Proposition [2] in Appendix [S] 
implies the following. 



lr The condition r > 3/2 ensures that functions within the kernel's native space are continuous on 
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Result 2: First-order discrete operators with finitely smooth kernels. Suppose cf> satisfies (26) with 
t > 3/2 + 1, define s = r - 1/2, and let /? satisfy s > /3 > 2. Then the following holds for all / G H^(M) 
and f € H^(M): 

||Gm/-V m /||l 2 (m) < Oih^p- 11 ), ||Gm/-Vm/|U oo( m)<0(^-V~^). 

n^Mf-VM-fiu^M) < o(/^-y-^), p M f-v M -f|| ioo(M) <o{h?- 2 p s -p). 

In addition, for very smooth target function^ these kernels give the following error rates: 
\\Guf - Vm/||l 3 (m) < C(/i 2s_1 ), || Gm/ - Vm/||l oc (m) < 0(h s ~ ), 

Uoc(M) 



|Aif-VM-f||wM) < 0{h 2s ^), \\D M t-V M -t\\ L , M) <0(h 2s - 2 ) 



Finding bounds on the discrete Laplacian is less straightforward - the repeated process of "interpolate, 
then differentiate" complicates matters. This will ultimately result in more factors of p in the estimates, 
even in the case of very smooth target functions. Unfortunately error bounds for the discrete Laplacian 
constructed from infinitely smooth kernels are currently unavailable. However, we will be able to provide 
high-order convergence rates in the discrete Laplacian for any kernel satisfying ( 26 1 . The following comes 
from Theorems [5] and [7J Appendix [A] 



Result 3: Discrete Laplacian with finitely smooth kernels. Suppose cj> satisfies (26) with r > 3/2 + 2, 
define s = r - 1/2, and let (3 satisfy s > j3 > 3. Then the following holds for all / G H^jM) and f G H^(M): 

||Lm/-A m /|U 2( m) < 0(h^ 2 p 2 ^ +1 ), ||Lm/-Am/||l oo( m) <0(^V (S ^ )+1 )- 

In addition, for very smooth target functions^ these kernels give the following error rates: 

||L M /- A m /||l 2 (m) < 0(h 2s - 2 P ), ||i M /-A M /|| ioo(M) <0(h 2s ~ 3 p). 

For solving time-dependent PDEs, one often is only concerned with approximating derivatives at the 
nodes since typically these are the only locations the numerical solution is known. If our point set X is 
always chosen in such a way that p is bounded by a fixed constant and / is any target function in a Sobolev 
space smooth enough to guarantee continuous second-order derivatives (that is, / G H^(M) with (3 > 3), 
Result 3 guarantees that 

\\L x fx - (Am/JMIi. < \\Luf - A M f\\ Loo( u) < 0(h^ 3 ) -> 0. 

Now that we have established spatial accuracy, we need to discuss the stability of numerical PDE solvers 
based on these discrete differential operators. 



6 Eigenvalue stability 

A necessary condition for stability of the MOL approach described in Section [4] is that the eigenvalues of the 



differentiation matrices Lx in ( 24 ) must be in the stability domain of the ODE solver used for advancing the 
system in time. We propose using either fully-implicit backward differentiation formulae (BDF) methods 
in the case where there is no forcing terms, and semi-implicit BDF methods (SBDF) in the case of forcing 
terms (see Sections [7] and [8] for details). This means that, at the very least, all eigenvalues must be in the 
left half plane. Unfortunately, the construction of the discrete approximation of Am does not guarantee that 
this property will hold for Lx- Fortunately however, we have found from numerous numerical experiments 
that provided the mesh norm h is small enough (i.e. the surface is well-discretized), then it appears that all 
of the eigenvalues do in fact lie in the left half plane. To illustrate these observations, in Figures [2]ja)-(f) the 
eigenvalues of L x are plotted for the surfaces and node sets depicted in Figure [T] Results are shown for both 
the IMQ kernel and the v — 7 Matern kernel, which is C 10 (IR 3 ). While not presented here, similar results 
were obtained for the Matern kernels with v < 7. The shape parameters, which are listed in the figure, have 
been chosen for this experiment so that the condition numbers of ([5| are around 10 10 . We find that even 
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Figure 2: Eigenvalues of the discrete Laplacian for the surfaces and node sets shown in Figure [T] x's 
correspond to the IMQ kernel with shape parameters chosen as (a) e = 2.8, (b) e = 4, (c) e = 5, (d) e = 2.8, 
(e) e = 2, and (f) e = 6.5, while o's correspond to the v = 7 Matern (MA) kernel with shape parameters 
chosen as (a) e = 8, (b) e = 12, (c) e = 15, (d) e — 8, (e) e — 6, and (f) e = 16. 



when the nodes are quite scattered (as in the case of the Dupin's cyclide and Bretzel2) the eigenvalues all 
lie in the left half plane in these examples. 

For node sets that are not sufficiently dense, we have found cases where some eigenvalues of the differ- 
entiation matrix shift over to the right half-plane. To illustrate this, we repeat the above experiment for 
the "bumpy sphere" with only N = 1035 nodes instead of TV = 5256. Figure |3][a) shows this reduced node 
set on the "bumpy sphere" , while part (b) of this figure shows the corresponding eigenvalues of the discrete 
Laplacian for both the v — 7 Matern and IMQ kernels. The shape parameters here were chosen so that 
the condition numbers of ^ are similar to those from the N = 5256 experiments. We can see from this 
latter figure that the majority of the eigenvalues fall in the left half plane for both kernels, but that there are 
now a few that fall in the right half plane. Fortunately, these "unstable" eigenvalues appear to correspond 
to high frequency eigenvectors as illustrated in part (c) of Figure [3] for the Matern kernel (a similar result 
holds for the IMQ kernel and the results are thus omitted). Unphysical eigenmodes of this type also appear 
in other contexts of kernel methods for PDEs 25,26,28 . For these applications, Fornberg and Lehto 28 



describe a global hyperviscosity stabilization procedure for selectively shifting the unstable eigenvalues to 
the left half plane without a reduction in accuracy. Given the similarity of the unstable eigenmodes in the 
current problem, it seems reasonable to expect this hyperviscosity procedure to also be applicable in the 
present context. In the numerical experiments of this study that follow in the next section, all eigenvalues 

2 See the discussion preceeding Proposition |5] in Appendix |a] for details. 

3 See Theorem [7] for details. 



11 




(a) (b) (c) 

Figure 3: (a) "Bumpy sphere" with N — 1035 nodes, (b) Eigenvalues of the discrete Laplacian corresponding 
to part (a), x's correspond to the IMQ kernel with shape parameter e — 2, while o's correspond to the 
v = 7 Matern (MA) kernel with shape parameter e = 6. (c) Eigenvector corresponding to the eigenvalue of 
the discrete Laplacian with the largest real part shown in part (b) using the v — 7 Matern kernel. 

of discrete Laplacian lie in the left half plane. We leave exploration of stabilization via hyperviscosity to a 
separate study. 



7 Numerical convergence results 

In this section we demonstrate the convergence of our numerical method for approximating the scalar diffusion 
equation with forcing (23 1 on two surfaces: the unit sphere S 2 and the torus shown in Figure [ljd) . In all 
the tests we use the standard fourth-order backward differentiation formulae (BDF4) method for advancing 
the approximate solution in time. The time-step is set to At = 10~ 4 to ensure that spatial errors dominate 
over temporal errors. 

Results for three kernels of varying smoothness are presented to illustrate the different convergence rates 
that are possible. The first two are the v = 4 and v = 6 kernels from the Matern family ([6]), which are C 4 (]R 3 ) 
and C 8 (1R 3 ), respectively. The third kernel is the infinitely smooth IMQ (I7|. For the Matern kernels, the 
analysis from Section [5] applies so we expect an algebraic rate of decay in the approximate solutions. While 
we don't have estimates on the approximation of the surface Laplacian for the IMQ kernel, the estimates 
from Result 1 in Section [5] on the discrete surface gradient and divergence suggest that this infinitely smooth 
kernel should give rates of approximation that are faster than any polynomial order (provided the underlying 
target is sufficiently smooth). We fix the shape parameters of the kernels in all the tests as follows: e = 4 
for v = 4 Matern, e = 8 for v = 8 Matern, and e = 3 for IMQ. These values have not been optimized and 
are chosen simply to illustrate the convergence rates of our method. 

Finally, we measure the errors in the numerical solutions using both the discrete two-norm and max- 
norm. To define the discrete two-norm we use an approximation of the continuous L2(M)-norm based on 
quadrature rules defined at the collocation nodes. We use the following abuse of notation to denote our 
discrete two-norm: 

ll/lk(M) := (X [/(x)]2dX ) V2 " (l>[/( x *)] 2 ) : = WfW^ 

where {wi}fL 1 are quadrature weights for the surface integral at the collocation nodes {x^}^ on the mani- 
fold. These weights are computed using a second-order accurate midpoint rule approximation on the surface. 
We use the standard definition for discrete max- norm and denote it with the standard notation of . 
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(c) Matern v = & (d) IMQ 

Figure 4: (a) Initial condition for the forced diffusion problem on the unit sphere, (b)-(c) Convergence 
results for the forced diffusion problem for three different radial kernels. These last three plots are displayed 
on a log-log scale with the horizontal axis being the square root of the number of collocation nodes N, which 
satisfies \/N ~ 1/h. Black dashed lines in (b) and (c) are the lines of best fit to the last 4 values and the 
slopes of the lines are given in the legends. 



7.1 Diffusion on the unit sphere with forcing 



For the first example, we use a test problem first presented in 11 . For this test, an artificial solution to (23) 
is specified (with 5 = 1) and the forcing function / is chosen so that this solution is maintained for all time. 
The solution is given by 

23 

u(t, x) = exp(— 5t) ^ exp (— 10 arccos(£ fc • x)) , (27) 

k=l 

where £ fc , k = 1, . . . , 23 are randomly placed points on the surface of the sphere. The term in the sum is 
a Gaussian centered at £ fe with the distance measured using the geodesic distance. The solution is clearly 
C°°(S 2 ) and its value at t = is displayed in Figure [i| a). In our results, the forcing function corresponding 
to (|27|) is computed analytically and evaluated implicitly in time. 



The test calls for comparing the errors in the approximate solution of (23 1 at time t = 0.2 at various 



spatial resolutions. For our spatial resolutions, we use minimum energy (ME) node sets on the unit sphere 
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discussed in Appendix [BJ| The different sizes of the node sets are N = 1024,1444,2025,3136,4096, and 
5041, which are all perfect squares. Figures |5]Jb)-(d) display the results of the convergence study for these 
node sets and the three respective radial kernels. Since the mesh norm for the ME nodes satisfies h ~ 1/y/N, 
we plot the relative errors vs. \/N. All three convergence plots are on a log-log scale. 

Figures ^b) and (c) illustrate the algebraic convergence rate of the approximate solution for the two 
Matern kernels. The black dashed line in these figures shows the lines of best fit to the data and the slope of 
the different lines is given in the legend. Since the solution is C°°(S 2 ), we expect the estimates from Result 
3 to apply to the approximations of the surface Laplacian in this test problem. Note that node sets we are 
using are quasi-uniform, so the mesh ratio that appears in Result 3 will have little or no effect here. For 
the v — 4 Matern kernel this result predicts the approximate surface Laplacian to converge like 0(h 5 ) in 
the ^2-norm and 0(h 4 ) for the ^-norm, while for the v — 6 Matern kernel these convergence rates should 
be 0(h 9 ) and 0(h s ). We see from Figures [^b) and (c), however, that the measured rates of convergence 



for the numerical solutions to the forced diffusion problem ( 23 1 are much higher than the predicted rates 
for the surface Laplacian. One reason for these higher rates may be that the estimates from Result 3 apply 
to the entire manifold, while we are only measuring the errors at the collocation nodes (which are the only 
locations where the approximate solution is known). We may therefore be benefitting from some form of 
super-convergence, which is known to occur at the nodes in certain periodic spline methods |62| . A more 
detailed study of this phenomenon is currently under way [35| . 

Figure |4^d) displays the results for the IMQ kernel. The results in this plot suggest that the convergence 
of the approximate solution is converging at a rate faster than any polynomial order. Further numerical 
investigations not presented here seem to indicate that the solution is converging at an exponential rate. 
This type of convergence has also been observed for collocation methods based on infinitely smooth kernels 



for the scalar transport and nonlinear shallow water wave equations on the surface of the sphere 25 26 



7.2 Diffusion on a torus with forcing 



For the second example, we again consider computing the solution of (23 1, but with the torus shown in 
Figure [TJd) and described in Appendix B.4 Like the previous test, we specify the solution and then choose 



the forcing function so that this solution is maintained for all time. The solution is given by 

u(t, x) = i exp(-5t)x (x 4 - 10xV + 5y 4 ) (x 2 + y 2 ~ 60z 2 ) , (28) 
8 

where x is a point on the Torus; see Figure [5|a) for a plot of u(0,x). We compute the forcing function 
analytically and evaluate it implicitly in the BDF4 scheme. A key ingredient to computing the forcing 
function is determining the surface Laplacian of it(i,x). For convenience, we provide the equation for this 
below: 

A m u(t, x) = - ~ exp(-5i) [x{x 4 ~ \§x 2 y 2 + 5y 4 ) 

(10248/? 4 - 34335p 3 + 41359p 2 - 21320p + 4000)] , (29) 

where p = \J x 2 + y 2 . 

As in the previous test, we compare the errors in the approximate solution of (23) at time t = 0.2 at 
various spatial resolutions. We also use ME node sets on the torus for the spatial discretization as discussed 
in Appendix |B~4l The sizes of the different node sets are N = 500,750,1000,2000,3000, and 4000. The 
results of the convergence study for these node sets and the three radial kernels are displayed in Figures 
[5|b)-(c). As in the previous test, the mesh norm for the ME nodes on the torus satisfies h ~ 1/y/N, so we 
again plot the relative errors vs. y/N. 

Figures (5^b) and (c) display the results for the v = 4 and v = 6 Matern kernels, respectively. Like the 
previous test, we clearly see the algebraic rates of convergence for these kernels, and the measured rates 
of convergence of the approximate solutions for both kernels are much higher than the predicted rates of 
convergence for the discrete Laplacian from Result 3. We again believe these higher rates may be related 
to our restriction of measuring the errors at the collocation nodes. Finally, we note that the measured 
convergence rates for the current test problem and the previous one are comparable, with only a slightly 
higher observed rate for the current one. 
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(c) Matern v = 6 (d) IMQ 



Figure 5: (a) Initial condition for the forced diffusion problem on a torus, (b)-(c) Convergence results for the 
forced diffusion problem for three different radial kernels. These last three plots are displayed on a log-log 
scale with the horizontal axis being the square root of the number of collocation nodes N, which satisfies 
y/~N ~ l/h. Black dashed lines in (b) and (c) are the lines of best fit to the last 3 values and the slopes of 
the lines are given in the legends. 

Figure [5jc) displays the results for the IMQ kernel. The results indicate the approximate solution is 
converging at a rate higher then any polynomial degree, which is entirely inline with the previous results. 
While not presented here, a further analysis of the numerical results indicate that the rate of convergence is 
again exponential for this kernel. 

8 Applications 

In the applications we consider systems of PDEs of the form ([I]). We use the third-order, semi-implicit, 
backward differentiation formulae (SBDF3) method discussed and analyzed in II] for advancing the semi- 
discrete systems in time. This scheme treats the diffusion terms implicitly and the (non-linear) reaction 
terms explicitly. We bootstrap the simulation with one step of the first order SBDF followed by one step of 
the second order SBDF (see [I] for details on these schemes). Upon initial LU decompositions of the linear 
systems that result from the implicit discretization, each time-step requires 0(N 2 ) operations. 

We only present results for the IMQ kernel using the same shape parameters as used in the stability tests 
(refer to the caption of Figure [2] for these values). However, similar results were obtained when using other 
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shape parameters and when using the Matern kernels. 



8.1 Turing patterns 

Since Turing's classical paper [64] that suggested how certain non-linear models of reaction and diffusion 
can lead to stable, heterogeneous pattern formations, there has been an explosion of research in reaction- 
diffusion-type models for various kinds of morphogenesis 44 51 53 . Most of the theoretical and numerical 
investigations of these models have been restricted to flat planar domains, but there is growing interest in 
studying these models on more general surfaces to understand how curvature may effect pattern formation. 
Additionally, there has been interest in the computer graphics community for using Turing models to generate 



interesting textures and patterns on surfaces 65 



We illustrate how our method can be applied to these types of pattern formation problems by applying 
it to the Turing system (a linearized Brusselator model) from [7 67: 



du 

at 

dv 

— =S v A M v + fiv ^ 1 + ~^~ uv \ + u h + T % v ) ■ 

"V" 

fv(u, V) 



—5 u A M u + au(l — nil 2 ) + v(l — r 2 u), 

fu(u,v) 
OLT\ 



(30) 



Here it and v are morphogens with u the "activator" and v is the "inhibitor" . If a = —7 then (it, v) — (0, 0) is 
a unique equilibrium point of this system. By changing the diffusivity rates of it and v an instability can form 
that leads to different pattern formations. The cubic coupling parameter t\ favors the formation of stripes, 
while the quadratic coupling parameter T2 favors the formation spots |7j. The spot pattern formations are 
more robust than stripes and take far less time to reach "steady-state" . The model ( 30 1 and similar models 



have been previously used to illustrate the applicability of some of the other numerical methods for PDEs 
on surfaces (8|[9j[TTJ[50}[57[[58}[65] . 

Figure [6] shows the results from our numerical solutions of (30) on various surfaces using parameters that 



lead to both spot and stripe patterns. Similar to the experiments on the surface of the sphere in 67 , we set 
the initial values of it and v to random values between —0.5 and 0.5 in a thin strip around the "equator" of 
each surface and u = v = elsewhere. Values for all parameters were also motivated from the values used 
67 and are listed in Table [l] A time-step of At 



0.01 was used in all the simulations and the numerical 
solutions were computed until steady-state was reached. The resulting stripe and spot patterns in Figure [6] 
are qualitatively similar to those obtained from other numcr ical methods (8|[Qj[TTJ[50j[5Sj[65] . 



Surface /Pattern 




S v 




a 


P 


7 






Red blood cell/spots 


4.5 


10- 


-a 


0.899 


-0.91 


-0.899 


0.02 


0.2 


Red blood cell/stripes 


2.1 


10" 


-3 


0.899 


-0.91 


-0.899 


3.5 





Bumpy sphere/spots 


4.5 


10" 


3 


0.899 


-0.91 


-0.899 


0.02 


0.2 


Bumpy sphere/stripes 


2.1 


10" 


-3 


0.899 


-0.91 


-0.899 


3.5 





Dupin's cyclide/spots 


4.5 


10" 


2 


0.899 


-0.91 


-0.899 


0.02 


0.2 


Dupin's cyclide/stripes 


1.89 


■ 10 


-2 


0.899 


-0.91 


-0.899 


3.5 





Bretzel2/spots 


2.1 


1Q- 


-3 


0.899 


-0.91 


-0.899 


0.02 


0.2 


Bretzel2/stripes 


8.87 


■ 10 


-4 


0.899 


-0.91 


-0.899 


3.5 






Table 1: Values of the parameters of (30) used in the numerical experiments shown in Figure [6] In all cases 
S u = 0M6S V . 



8.2 Spiral waves in excitable media 

Spiral waves can be observed in many excitable chemical, biological, and physical media |18||43 | [66) . Impor- 
tant examples include Belousov-Zhabotinsky chemical reactions and electrical activity in the membranes of 
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Figure 6: Turing spot and stripe patterns computed from the model ( 30 ) on various surfaces. The pseudocolor 
plots are for the activator u once steady state is reached. In all plots red corresponds to a high concentration 
of u and blue to a low concentration. For all simulations the time-step was set to At = 0.01; parameters 
used to get the different patterns are given in Table [I] 



organisms. While numerical methods have been developed for these models in planar domains (see [5j[6] and 
the corresponding software EZ-Spiral), there has been growing interest in studying these models on non- 
planar surfaces since most physically relevant problems occur on curved surfaces and curvatrue can effect 
[BJ[40](52] 



the wave processes 



To illustrate how our kernel method may be applied to these models, we focus 
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Figure 7: Snap shots of spiral wave patterns computed from the Fitzhugh-Nagumo model (31) at different 



times. The left column shows the excitation variable u on the bumpy sphere at times t = 42.24,42.80,43.76 
and 44.8. The right column shows u on Dupin's cyclide at times t = 42.76,43.48,44.40 and 45.24. The node 
sets used in the simulations are shown in Figure [I] (c) and (e) , respectively, and the time-step was set to 
At = 0.02. The colors in the pseudocolor plots change linearly from white (u = 0) to black (it = 1).. 

on a simple two- variable reaction-diffusion model used in |5j: 

du 1 / v + b 

— =d u A u u + —u (1 - u) \ U 

at a \ a 



fu(u,v) 



dv 

di 



(31) 



=(5„Amv + u — v , 

fv(u,v) 
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where u and v can be viewed as some chemical concentrations or as membrane potential and current. The 
parameters a, b, and a govern the reaction kinetics and 5 U and 5 V are the diffusivities of the u and v species 
respectively. The parameter a is chosen as a << 1 so that the u field takes on the values u = or u = 1 



almost everywhere, with a thin interface (or reaction zone) separating these two regions. The system (31) 



is of Fitzhugh-Nagumo type [22] and provides a simple model for dynamics of many excitable media. These 



types of systems have also been studied on the surface of the sphere 38 72 . 

Figure [7] displays snap shots from the numerical solution of (31 ) on the bumpy sphere and Dupin's cyclide 



at different times of the simulation. For all of these experiments we set the initial conditions to 

tt(0, x) = - [1 + tanh(2x + y)] , 
u(0,x) = - [1 - tanh(3z)] , 

and At = 0.02. Values for the parameters that remain fixed for the two surfaces are as follows: a — 0.75, 
b = 0.02, a = 0.02, and S v — 0. The values of S u vary depending on the surface as follows: 5 U — 1.5(27r/50) 2 
for the bumpy sphere and S u = 2.5(27r/50) 2 . These values are similar to those used in one of the experiments 
from [6] on a 2-D planar domain. The two quasi-periodic, counter-rotating spiral waves seen on both surfaces 
in Figure [7] are qualitatively similar to those observed on the sphere in [38] . These spiral waves remain intact 
but meander around the surfaces throughout the simulation time. 

9 Concluding remarks 

In this paper, we have introduced a new kernel method based on collocation with RBFs for constructing 
discrete approximations to differential operators on surfaces, with a focus on surfaces of dimension 2 em- 
bedded in R 3 . The method is different than other methods for approximating surface differential operators 
in that it only requires "scattered" nodes on the surface and the corresponding normal vectors, it does not 
require expanding into the R 3 embedding space, and it can give high-orders of accuracy. We have studied 
the approximation power of the method and have provided error estimates for target functions of various 
smoothness. We have used the method to construct discrete approximations to the surface Laplacian and 
combined it in a method-of-lines approach to numerically solve diffusion and reaction-diffusion equations on 
surfaces. We presented numerical results for two forced diffusion equations, illustrating the convergence of 
the method for kernels with finite and infinite smoothness. The numerical results for the finitely smooth 
kernels indicate that the convergence rates are faster than our theory predicts, which may be related to 
some form of super-convergence common to spline methods. Finally, we have illustrated the flexibility of our 
method by successfully applying it to two important problems in biology and chemistry. 

The primary shortcoming of the method is the 0(N 2 ) computational cost for applying these operators. In 
a follow up study, we intend to address this issue by using a local RBF finite-difference (RBF-FD) method 
for constructing the approximate surface differential operators. This method has been used successfully 
in 24 28 to overcome the computational cost of the global method 25 26 for the shallow water wave 
equations on the sphere and has recently been adapted to multi-CPU architectures in 10 . We expect the 
cost of constructing and applying these discrete operators based on RBF-FD to be reduced to O(N), with a 
minor decrease in the order of accuracy. These enhancements may make the method competitive for solving 
diffusion and reaction-diffusion equations on evolving surfaces, which is also important in many biological 
applications. 

A Convergence Results 

In this section we present the convergence results for our discrete differential operators. The results presented 
here are for the special case when M C R 3 is a two dimensional surface, although similar results hold in 
higher dimensions. The arguments will depend h eavi ly on error estimates for kernel interpolants on smooth, 



embedded submanifolds of R given recently in 36 . We will frequently reference results from that paper 



throughout this section. Sobolev error estimates for the surface gradient and divergence operators will 
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immediately follow the results in 36 . However, the estimates for the discrete surface Laplacian will require 
more work. 



Proposition 1. Let 1 < q < oo, (f> be a positive definite kernel satisfying (25) with r > 3/2 + 1, and define 
s = t — 1/2. Then there is a constant hf& depending only on M such that if a finite node set X C M satisfies 
h<h u , for all f e 7V^(M) and f e (A^(M)) 3 we have 

||Gm/-v m /||l 9 (m) < CTi 8 - 1 - 2 ^ 2 - 1 /^!!/!!^), 

||£>Mf-V M -f|U a (M) < C/ 1 s - 1 - 2 < 1 /2-l/9)+ ||f||^ (M) , 
where (x) + = x if x > and is zero otherwise. 
Proof. We will focus on the discrete gradient operator. We have 

||Gm/ — Vm/||l,(m) = ||Vm^/ — Vm/||l,(m) < C\\I$f — /||vki(M)) 
where W^M) is the L q Sobolev space of order 1 (see Section 2, |36|). The assumptions on r allow us to 



apply Theorem 4.6 in 36 , and the results follow. The error bounds in Theorem 4.6 in [36] carry over to the 



vector-valued case, thus bounds for the divergence operator error are obtained in a similar manner. □ 



When 4> satisfies (26), i.e. cj) has finite smoothness, the estimates come in two types: the first concerns 
targets that may be too rough to be within the native space, and the second applies to functions with 
additional smoothness. This additional smoothness is measured with the inverse of the integral operatoij^] 

Tf(x) := f <p(y,x)f(y)d[i(y). 
Jm 

We denote a vector version of this operator by T, which simply applies T to each component of a 3- 
dimensional vector field. Similar to the proof above, the result below follows from Theorem 4.12 and 
Corollary 4.10 in [36]. 



Proposition 2. Let 1 < q < oo, and let <f> be a positive definite kernel satisfying (26) with t > 3/2 + 1, and 
define s = r — 1/2. Then there is a constant Iim depending only on M such that if a finite node set X C M 
satisfies h < hja, we have the following: 



1. Rough target functions. Let f3 be such that s > /3 > 2. Then for all f e iF(M) and f € H' 3 (M) 
we have 

I|Gm/-v m /||l 9 (m) < c^- 1 - a(1/2 - 1/9 V _/ 'll/IU/»(M)> 

||Awf-V M -f|U 9(M) < C/^-^-^V^IIfllH^M). 

2. Smooth target functions. For all f € L 2 (M) such that T~ l f e L 2 (M) and f € L 2 (M) such that 
T~He L 2 (M), we have 

||Gm/-V m /||l 9 (m) < ChP-i-W-W+WT-^WivW, 
||-D M f-V M -f|| Lg(M) < CT/» a — 1 - 3 < 1 / 2 - 1 /«)+ ||T- X f Hx.^. 

Obtaining convergence rates for the discrete Laplace operator is not as straightforward, and will require 
extra tools. The first is the "zeros lemma" from Narcowich, Ward and Wendland [54 . We will use the 
manifold version stated in [36] Lemma 4.5], with parameters adapted to our situation. 

Proposition 3. Let 1 < q < oo. t e K with t > 1. Let fi e N satisfy < \x < [t - 2(1/2 - l/q)+] - 1. Also, 
let X C M be a discrete set with mesh norm hx,M- Then there is a constant depending only on M such that 
if hx,m < Cm and ifuG H (M) satisfies u\x = 0, then 



4 With this integral o pera tor comes the pseudodiflerential operators T r , r > 0. A function / is in the native space if and 
only if T~ 1 ' 2 / £ L%(M.) 36| Proposition 4.9], so we expect functions such that T -1 / £ Z/2(M) to be twice as smooth. 
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Proposition [3] also holds with full Sobolev norms, and a similar result holds for vector fields. 

Since the operator Lm is a mixture of differential and interpolation operators, it will be beneficial to 
bound the norm of the interpolation operator and its associated error in several different Sobolev spaces. In 
particular, we will use the following. 



Lemma 4. Let <j) satisfy {26) with r > 3/2, and define s = r — 1/2. Let be such that s > j3 > 1. Then 
there is a constant h$/i depending only on M such that if a finite node set X C M satisfies h < hwi then for 
all f e H p (M) we have 

IIW - /llff£(M) ^ Cp S_/3 ll/ll Hf>(M) an d II VII-H^M) < Cp S ~^ll/l|jf3(M)- 
Proof. The last estimate in the proof of Theorem 4.12 in 136] is 

11/ - WIIh^cm) ^ Cp^^ll/IU^CM)- 

Applying a triangle inequality and observing that 1 < p, we get 

ll^/ll_ff 3 (M) < 11/ - ^/IIh^(m) + II/IIh^(m) < c(p s_/3 + i)ll/llif 3 (M) < Cp s_/3 II/IIh^(m)- 



□ 



Now we are poised to present the following theorem. 



Theorem 5. Let 1 < q < oo, <f> satisfy (26) with r > 3/2 + 2, and define s = r — 1/2. Let ft be such that 
s > p > 3. TTien t/iere is a constant hyi depending only on M such that if a finite node set X C M satisfies 
h < fei, then for all f £ H /3 (M) we have the estimate 

\\Lmf ~ Am/|U, (m) < Ch^-W 2 - 1 '^ p^ s -V +1 \\f\\ Hmy (32) 

Proof. First, we have 

||£m/ — Am/|U,(m) < II^m/ — Am^/|U,(m) + II Am/ — Am-^/Hl^m)- (33) 
For the rightmost term, the assumption /3 > 3 allows us to use [36] Theorem 4.12] to get: 

||A m /-AmVI|l,(m) < C\\f-I^f\\ wsm < C^-^-^+p^H/Uff^M). 
For the other term in ( |33[ ) , we have 

II-^m/ — Am^/Hl^m) = ||Vm • (^$(Vm^/)) — Am/^/Hljm) 

= ||Vm • (Lj>(VuL^f)) — Vm ■ (Vm^/)||l,(m) 

< C||/$(Vm^/) — Vm^/||wi(m)- (34) 

Note that the last estimate is the interpolation error for the target function g := Vm^/- 

We claim that g is in H*(M) for all t < 2s - 2. Since satisfies ([26]), I^f G i^(K 3 ) for all v < 2r - 3/2. 
By the trace theorem restricting to M puts I^f in i?"- 1 / 2 (M). Thus g G IF~ 3 / 2 (M) for all < 2r - 3/2, 
which is equivalent to g G H*(M) for all t < 2s — 2. In particular, g G H^ _1 (M). Also, since (3 > 2 then we 
have 

1 < \{fi - 1) - 2(1/2 - 1/?)+] - 1, 

and that j3 — 1 > 1. This allows us to estimate p4| with Proposition [3] (with parameter /z. = 1 and target 
smoothness t = f3 — 1) to get 

||I*(VmW)-VmWIIwj(m) < cM' 3 - 1 )- 1 - 2 ( 1 / 2 - 1 ^+||/$g-g||H^i(M)- 

= C7^- 2 - a (Va-i/«) + ||j, g _ g || Hfl _ 1(M) . 

Further, since j3 — 1 > 1, we can also apply Lemma [4] and with two applications of it we get 

||J*g-g||H0-i(M) < Cp S "^ +1 ||g||H-3-i(M) = C"p S "^ +1 1| V M /0/|| H /3-i(M) (35) 

< Cp^ +1 ||^/||^ (M) <c P 2 ^ +1 \\f\\ H , m . 

This completes the proof. □ 
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Continuing with our error analysis, we now shift our attention to target functions that are very smooth. 
First we will need a lemma. 



Lemma 6. Let <j> satisfy (26) with r > 3/2, and let s = r — 1/2. Then for all f G L 2 (M) such that 
T~ 1 f G L 2 (M) we have the following estimate 

11/ - Wllif=(M) - Ch s \\T~ 1 f\\ L2 ( M - j . 

Also, for all f G L 2 (M) such that T _1 f G L 2 (M) we teue 

||f-/*f||H«(M) < C^HT-^Hl^m). 



Proof. The first estimate is established in the proof of Corollary 4.10 in 36 . By working component-wise 
the second estimate follows from the first. □ 



Theorem 7. Let 1 < q < oo, f> satisfy (26) with r > 3/2 + 2, and define s = r - 1/2. Lei / G L 2 (M) fee 
smc/j that T~ 1 f G L 2 (M) and T _1 Vn/ G L 2 (M). Then there is a constant /ijj depending only on M suc/i 
i/iai i/ a finite node set XcM satisfies h < hj&> we have the estimate 

\\L M f - A M f\\ Lq{M) < Ch 2s - 2 - 2 W-Vi) +p (\\T-if\\ L2(u) + HT^VM/lkcM))- 

Proof. We begin as in the proof of the previous theorem. We have 

\\Lwif — Am/|U„(m) < H-^m/ ~ ^-M.I(j>f\\L q (u) + ||Am/ — A M /0/|| L(! ( M ). (36) 



To estimate the rightmost term we can use the error estimates in 36 Corollary 4.10] to get: 

||A m /-AmWHl,(m) < C\\f- Wlk|(M) <Ch 2s - 2 -^ x ' 2 - x l^+\\T- x f\\ L2{m) . 

For the other term in p6| , we proceed as in the proof of Theorem [5] to get 

\\Luf — Am^/|U s (m) < C||i$(g) — g||wj(M)> (37) 

where g = Vm^/, and as before we know that this function is in H*(M) for all t < 2s — 2. In particular, 
g G H S_1 (M) and we can estimate (37) with Proposition [3] (with t = s — 1 and fj, = 1) to get 



|/*(g) " gllwJCM) < C/ l S - 2 - 2 ( 1 /2-l/9) + || /<& g_g|| Hs 



L (M)- 



This is where the proof detours from that of the previous theorem. To estimate ||7$g — g||H s - 1 (M)i note that 
it is bounded it by following quantity: 

||7*(Vm^/) - 7*(Vm/)||h=- 1 (m) + ||7#(Vm/) — Vm/Hh^-^m) + ||Vm/ — Vm^/||h s -!(m) ■ 



II III 



We will bound each term individually. First we concentrate on I. An application of Lemma |4j which we 
may apply since s — 1 > 1, gives us 

||-f$(VM-f<?>/) - ^$(Vm/)||h s - 1 (m) — ||^$(Vm^/ — Vm/)!!^-^!^) 

< Cp\\S7mI<j,f - Vm/Hh^-hm) = Cp(III). 

Thus a bound for I LI will result in a bound for J. To bound LLL, we can apply Lemma [6] to get 

HVm-W - Vm/|h=-i(m) < cil-W - /||h s (m) < /j. s ||r~ 1 /|| J c 2 (M)- 

To bound //, we again employ Lemma [6j 

||7*(Vm/) — Vm/||h«-i(m) < ||^*(Vm/) — V m /||h s (m) < C^ s ||T _1 Vm/||l 2 (m)- 
This completes the proof. □ 
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B Surfaces and node sets from the numerical experiments 



B.l Unit sphere 

This manifold is, of course, described implicitly by 

M = {x = 0, y, z) G M 3 \x 2 + y 2 + z 2 = 1 } . (38) 
The node sets we use for discretizing the unit sphere are the minimum energy (ME) node sets of Womersley 



and Sloan 69 . These node sets are approximately uniformly distributed over the surface of the sphere and 



have the nice property that the mesh norm h and the separation radius q decrease uniformly like the inverse 
of the square root of the number of nodes N, i.e. h,q ~ Additionally, the nodes in these sets are not 
oriented along any vertices or lines, which emphasizes the ability of our method to handle arbitrary node 
layouts. They have been used quite successfully in many other RBF applications, e.g. 25,26,28,33,34,55,70 . 



B.2 Red blood cell 

This manifold is a mathematical model for human red blood cells in static equilibrium conditions. The 
model was first proposed in [l9j and has been used in many subsequent studies (e.g. [48])- The model can 
be described parametrically as follows: 



(x,y,z) e 



r cosAcos£?, y = r sin Acos#, z 



1 . 

- sin I 
2 



(c + c 2 cos" 6 + C4 cos 



(39) 



where -tt/2 < 9 < it/2, -tt < A < tt, r = 3.91/3.39, c = 0.81/3.39, c 2 = 7.83/3.39, and c 4 = -4.39/3.39. 
The node sets we used for discretizing this manifold were obtained using a radial projection of the ME points 
for the surface of the sphere described above. 



B.3 Bumpy sphere 

This manifold is constructed from the "bumpy sphere" surface [3] using the following procedure. First, 
N = 5256 points on the bumpy sphere surface were obtained from http://shapes.aimatshape.net. This 
surface is homeomorphic to the unit sphere, so spherical coordinates for each of the N = 5256 points on 
the bumpy sphere were obtained. Upon normalizing the original nodes by the maximum distance of any 
of the nodes from the origin, a parametric model for the surface was constructed using the RBF geometric 
modeling technique from [61] . The original (normalized) N — 5256 points on the bumpy sphere are used as 
the computational nodes in all the numerical experiments and the parametric model is used for computing 
the normal vectors to the surface. 



B.4 Torus 

This manifold is given by the implicit equation: 

= (x, y, z) G M 3 



(i- V* 2 + y 2 ) 2 + 



(40) 



The node sets used for discretizing this manifold were obtained by arranging the nodes so that their Reisz 
energy (with a power of 2) is near minimal as described in Hardin and Saff 's seminal article 41 . As with 



the ME sphere nodes, these ME torus nodes sets provide a near uniform discretization of the manifold with 
an approximately uniform decrease in the mesh norm h and the separation radius q. The node sets used in 
the experiments were kindly given to us by Drs. Douglas Hardin and Edward Saff and Ms. Ayla Gafni from 
Vanderbilt University. 
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B.5 Dupin's cyclide 

This manifold is given by the implicit equation: 

M= {x= (x,y,z) €E 3 (x 2 + y 2 + z 2 - d 2 + b 2 f -A{ax + cdf - 4b 2 y 2 = 0} , (41) 

where a — 2, b — 1.9, d = 1, and c 2 = a 2 — 6 2 . The node set for discretizing this manifold was obtained 
from the algorithm of Palais, Palais, and Karcher (PPK) [56], which will be included in a future version 
of the amazing mathematical visualization software package 3-D-XplorMath lj. This algorithm generates 
"uniformly random" point clouds for surfaces and works for both implicit and parametrically defined surfaces. 
To obtain the node set displayed in Figure [TJe) we used the PPK algorithm to generate a point cloud 
consisting of 9562 points. We then thinned this point set by removing points that were too close together 
in order to increase the separation radius. The final node set ended up at N = 4948. Both choices for the 
starting and ending number of nodes were chosen somewhat arbitrarily. 

B.6 Bretzel2 

This manifold is given by the implicit equation: 



jx = (x, y, z) G E 3 (x 2 (l - x 2 ) - y 2 ) 2 + \z 2 = -L J 



(42) 



Unlike the the other surfaces, this one does not have an explicit parameterization. The node set for this 
domain was also obtained from the PPK algorithm. In this case we started with a point cloud of 10278 
points and thinned it to end up with a node set of N = 5041 nodes. 
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